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Abstract —In the present work, we consider a prototypical 
example of a PT-symmetric Dirac model. We discuss the 
underlying linear limit of the model and identify the threshold 
of the VT -phase transition in an analytical form. We then focus 
on the examination of the nonlinear model. We consider the 
continuation in the PT-symmetric model of the solutions of the 
corresponding Hamiltonian model and find that the solutions can 
be continued robustly as stable ones all the way up to the PT- 
transition threshold. In the latter, they degenerate into linear 
waves. We also examine the dynamics of the model. Given the 
stability of the waveforms in the PT-exact phase we consider 
them as initial conditions for parameters outside of that phase. 
We find that both oscillatory dynamics and exponential growth 
may arise, depending on the size of the corresponding “quench”. 
The former can be characterized by an interesting form of bi¬ 
frequency solutions that have been predicted on the basis of the 
SU(1,1) symmetry. Finally, we explore some special, analytically 
tractable, but not PT-symmetric solutions in the massless limit 
of the model. 

Index Terms —Nonlinear dynamical systems, nonlinear differ¬ 
ential equations, bifurcation. 

I. Introduction 

The study of open systems bearing gain and loss (especially 
so in a balanced form) is a topic that has emerged over the 
past two decades as a significant theme of study 0-0- 
While the realm of VT -symmetry introduced by Bender and 
collaborators was originally intended as an alternative to the 
standard Hermitian quantum mechanics, its most canonical 
realizations (beyond the considerable mathematical analysis 
of the theme in its own right at the level of operators and 
spectral theory in mathematical physics) emerged elsewhere 
in physics. More specifically, in optical systems 0.0 the 
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analogy of the paraxial approximation of Maxwell’s equations 
and of the Schrodinger equation formed the basis on which 
the possibility of VT -symmetric realizations initially in optical 
waveguide experiments was proposed and then experimentally 
implemented (6). The success of this program motivated 
further additional initiatives in other directions of experimental 
interest, including, but not limited to, VT- symmetric elec¬ 
tronic circuits 0 0 mechanical systems (9]] and whispering- 
gallery microcavities ]T0| . 

Another theme of research that has been receiving in¬ 
creasing attention recently, both in the physics and in the 
mathematics community is that of the nonlinear Dirac equa¬ 
tions (NLDEs). While such models were proposed in the 
context of high-energy physics over 50 years ago CD’ GT 
they have, arguably, been far less widespread than their non- 
relativistic counterpart CD the nonlinear Schrodinger (NLS) 
equation 0 > 0 . In recent years, however, there has been a 
surge of activity around NLDE models fueled to some extent 
by analytical solutions and computational issues arising in 
associated numerical simulations 0-0’ as well as by the 
considerable progress achieved by rigorous techniques towards 
aspects of the spectral, orbital and asymptotic stability of 
solitary wave solutions of such models p0| -p5] and towards 
criteria for their spectral stability (26), | |27| Although our em¬ 
phasis herein will be on the so-called Gross-Neveu model [ |28| 
(sometimes also referred to as the Soler model |[29]), we 
also mention in passing that another main stream of activity 
in this direction has been towards the derivation of NLDEs 
in the context e.g. of bosonic evolution (30) , CD (or light 
propagation (32) , (33) ) in honeycomb optical lattices. In the 
latter context, contrary to what we will be focusing on below, 
there is no nonlinear interaction between the fields (of the 
two-component spinor). 

Our aim in the present work is to connect these two 
budding areas of research, namely to propose a prototypical 
PT-symmetric nonlinear Dirac equation model (PT-NLDE). 
Motivated by the considerable volume of activity, as well 
as analytical availability of solutions within the Hamiltonian 
limit, we will focus on the Gross-Neveu (or Soler) model. In 
the next section, we will present the mathematical formulation 
of the model. We will analyze its linear limit and discuss the 
existence of a PT-symmetry breaking critical point, i.e., the 
point of a VT- phase transition. Then, we will turn to the 
nonlinear variant of the model exploring the conditions for the 
existence of a standing wave solution, as well as discussing 
the linear (spectral) stability setup for such a solution. Finally, 
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we will briefly touch upon the fate of the conservation laws, 
such as the power (squared L 2 norm) and the energy. As a 
remarkable feature, we find that the energy remains invariant 
within the nonlinear PT-NLDE model, a feature that certainly 
distinguishes the model from its NLS counterpart. In Section 
III we will examine the numerical properties of the standing 
wave solutions. Remarkably, we will find that these standing 
wave solutions are stable throughout their interval of existence 
tending to a linear limit of vanishing amplitude as the linear 
threshold (i.e., the threshold of the underlying linear model) 
of the PT-transition is approached. Since no instability is 
encountered within the exact PT-phase, we consider the 
propagation beyond the relevant critical point (i.e., under a 
quench in the gain-loss parameter 7 ) finding (a) the possibility 
of oscillatory motion that we identify with a bi-frequency 
state and connect to the invariances of the model and (b) 
the possibility of exponential growth. Lastly, a special case 
of vanishing mass solutions is analytically identified and also 
numerically examined. The remarkable feature in this case 
is that these solutions do not respect the PT-symmetry. In 
Section IV, we summarize our findings and present some 
interesting directions for future study. 

II. Model equations 

The system of choice in the present context will be the 
Gross-Neveu model in its generalized PT-symmetric (PT- 
NLDE) form: 

i d t U = d x V - g(\U\ 2 - \V\ 2 ) k U + mU + i 7 V, (la) 
idtV = -d x U + g{\U\ 2 -\V\ 2 ) k V - mV + v)U. (lb) 

We will restrict our considerations to the one-dimensional 
setting, as is evident from the above. Eqs. 0 are PT- 
symmetric because they are invariant under the transformation 
P: x —y — x , U —y U, V —y —V and P: t —y —t , i —y —i, 

U —y U, and V V. This transformation assumes that 
U(x,t) is spatially even and that V(x,t) is spatially odd. 
Without loss of generality, we will set the mass m = 1, 
except when explicitly indicated otherwise, while we will 
also use g = 1 (the coefficient of the nonlinear term can be 
scaled out, hence, when it is present, we only care about its 
sign). The key addition in this model in comparison to the 
earlier proposal of fl 8 | is the inclusion of the gain-loss term 
proportional to 7 in the (implicit) form of the Dirac matrix 75 
(cf. [ 11) multipliying the spinor ([/, V) T . In our case of two- 
component spinors, the role of 75 is played by the Pauli matrix 
07 . We note in passing that in its linear form, the model can 
be converted under a suitable transformation (associated with 
the so-call ed C-oper ator) to a Hamiltonian one with a reduced 
mass of \Jrn 2 — 7 2 [ 1 ], |34J. Although the linear version of 
the model was proposed and analyzed in these works, to the 
best of our knowledge, there has not been any previous work 
in the realm of the nonlinear variant i.e., of the PT-NLDE. 

It is straightforward to see that in the linear case (of g = 0), 
plane waves U = Ae l ( KX - ut ) and V = are solu¬ 

tions provided the dispersion relation uo = ±^m 2 + n 2 — 7 2 
is satisfied. Not only does the above formula have the 
characteristic Dirac form, but it also is consistent with the 


equivalence of the linear PT-Dirac equation with effective 
mass m = \Jm? — 7 2 , as per the above discussion. 

Having examined the linear states (plane waves) of the 
model, let us now turn to the nonlinear ones, and more 
specifically to standing waves. The relevant coherent structures 
will be of the form: 


U(x,t) = exp(—iAt)u(x), V(x,t) = exp(— iAt)v{x ). 

( 2 ) 

Once such standing wave solutions are calculated (e.g., by 
fixed point methods as will be discussed in the next section), 
their linear stability is considered by means of a Bogoliubov- 
de Gennes linearized stability analysis. We note here, in 
passing, that unfortunately, contrary to what is the case for 
the Hamiltonian limit of the model with 7 = 0, for which 
explicit solutions exist as: 


u(x)=\ 


v(x)=sgn(x) 


(m + A) cosh 2 (k/3x) 

(k +1)/3 2 

l/k 

Y m + A cosh( 2 k/3x) 

g 2 {m + Acosh (2 k/3x)) 



(3a) 


(m — A) sinh 2 (k/3x) 

(k + 1 )/? 2 

Y m + A cosh( 2 kf3x) 

_g 2 {m + Acosh (2 k/3x)) _ 


(3b) 


(where /? = \Jm 2 — A 2 ) in the present PT-NLDE we have 
been unable to identify such explicit solutions (with a notable 
exception for m = 0 discussed separately below). I.e, with 
the exception of the mm 0 solutions given below, we have 
not been able to identify additional analytical solitary wave 
solutions of the case with 7 7 ^ 0. In the same vein, it does not 
appear to be straightforward to generalize the transformation 
of (lj, p4| to the present nonlinear setting. 

We now consider the linearization of the standing wave 
solutions in order to extract information about their spectral 
stability. More specifically, considering small perturbations [of 
order 0(5), with 0 < 5 1] of the standing wave solutions, 

we substitute the ansatz 


U(x,t) = e lAt u(x) + 5(ai(x)e luJt 
V(x,t) = e~ iAt Ux) + S(a 2 (x)e [ujt 


hixye-^) jC4a) 
b 2 {x)*e-™* t ) l(4b) 


into Eqs. {!]), and then solve the ensuing [to 0(5)] eigenvalue 
problem cc(ai, h \, = M{a\, b\, b%) T with M 

being 


( L x L 2 \ ( J 0 \ 

M=\ - i 7 (5) 

V -£2 -LI ) V 0 J J 
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and 




g(\u \ 2 — |^| 2 ) fe — m + A 


- 0, 


d x — g(\u \ 2 — \v\ 2 ) k + m + A 


2 i ,|2\/e-i 


+gk {\uY-\vY) 


\u\ —uv 


-u V \v\ 


2 i ,\2\k-t 


L 2 =gk{\uY-\v\ z ) 


u —uv 


-uv V* 



where / is the identity matrix. The potential existence of 
an eigenvalue with non-vanishing real part (i.e., an eigenfre- 
quency c j with non-vanishing imaginary part) suggests the 
existence of a dynamical instability. If all the eigenvalues 
are imaginary (i.e., all eigenfrequencies cv are real), then the 
solution is spectrally (neutrally) stable. 

Finally, once the exact waveforms and their linear stability 
are identified, the corresponding full nonlinear dynamics of the 
scheme is monitored by means of the solution of Eqs. 0 in 
our numerical considerations of the next section. A natural, 
theoretically motivated aspect to consider in that regard is 
the (potential) preservation, by the numerical scheme, of the 
different conservation laws. To that effect, we examine the fate 
of the prototypical conservation laws (such as the power and 
the energy) in the context of our PT-NLDE model. 

Based on the power density, 

p(x) = \U(x,t)\ 2 + \V(x,t)\ 2 , ( 8 ) 

we can define the charge (or power, as it is also referred e.g. 
in the context of optics) 



From the dynamical equations 0 it is straightforward to 
show that the charge is not preserved. Instead, the following 
“moment equation” is satisfied: 


dQ 

dt 


4 7 J Rc(UV*)dx. 


( 10 ) 


Note that in the case of a standing wave state, dQ/dt = 0 and 
charge is conserved. 

Although the charge is not generally conserved, remarkably 
there is a conserved quantity in the form of the energy: 


E = 


= J [Re(C/* 


d x V-V*d x U) + m(\U\ 2 -\V\ 2 ) 


■kTi^r-w 


21 fc+1 


dx . 


(ID 


Notice that there is no 7 dependence in this formula. That 
is, this is the same definition for the energy as in the 7 = 0 
limit. But, intriguingly, dE/dt = 0 even for 7 / 0. This 


is rather unusual in our experience in PT-symmetric models 

) and is effectively related to the special form of introducing 
VT- symmetry through the 75 matrix. We note that in this 
form, it is not transparent (as it is e.g. in Schrodinger VT- 
symmetric models 11 ]) which component corresponds to the 
gain and which one to the loss. Effectively, isolating the time- 
dependence and the 7 -dependent term in the equations (i.e., 
id t U = i'yV and id t V = ijU and momentarily ignoring 
rest of the terms), it appears as if both components bear 
both gain and loss. Nevertheless, this 7 -independent effective 
energy conservation that we will numerically corroborate 
below is certainly worth additional examination in order to 
(7determine its origin. 

III. Numerical results 
A. Massive NLD equation with k = 1 

In the numerical computations presented herein, we have 
utilized spectral collocation methods in order to approximate 
the spatial derivatives of Eq. 0 . As discussed in j35), the 
Chebyshev collocation is, arguably, the most suitable method 
for approximating the relevant derivatives as it gives a better 
spectral accuracy. However, because of the non-Hermiticity of 
the system, the implementation of fixed-point methods, such 
as the Newton method used herein, requires a high amount of 
computer memory and more than ^ 500 grid points which, 
in turn, poses implementation challenges. Furthermore, the 
method has the drawback that the double humped solitons 
(which occur for A < 1/3 when k = 1), cannot be well 
resolved (i.e. the humps cannot be observed) because of the 
Chebyshev collocation including more points at the edge of 
the system in comparison to the center. Consequently, in what 
follows, a Fourier collocation scheme has been implemented. 
In suitable limit cases, we have checked that the results are 
similar for the different implementations and that no extra 
spurious eigenmodes arise in comparison to the Hermitian 
case. We note that hereafter we will focus on the case of 
k = 1 for our numerical implementation. 

The first numerical result found by studying the standing 
wave solution is the VT transition point. We have checked 
that this transition takes place when 7 = 7 pt with 7 pt = 
s/m 2 — A 2 (as we have taken m = 1 , 7 pt = a/ 1 — A 2 ); 
notice that this is consonant with our analytical prediction 
from the previous section in the case of wavenumber u = 0 . 
Figure [T] shows the profile of a typical soliton with nonzero 
7 . Importantly, we have observed (see Fig. [2]) that the non- 
Hermiticity does not introduce instabilities to our system. 
The relevant spectrum features a zero eigenvalue of alge¬ 
braic multiplicity four and geometric multiplicity two. This is 
present in the spectrum of the linearized equation due to the 
U(l) symmetry and due to the translational symmetry, which 
are both preserved when 7 / 0 , hence both the algebraic 
and geometric multiplicity of this eigenvalue are preserved 
for all values of 7 , as is the presence of two generalized 
eigenvectors. The spectrum also features the eigenfrequencies 
uj = ±2A fl9) which can not leave the imaginary axis 
since their presence in the spectrum is due to the SU(1,1) 
invariance (see the more detailed discussion below), which is 







4 



X X 

Fig. 1. Real and imaginary (top and bottom panels) part of the spinor 
components (left and right panels) for A = 0.8 and 7 = 0.3. 



Fig. 2. Real part of the eigenfrequencies dependence for A = 0.3 (left) and 
A = 0.8 (right). Notice the existence on the left panel of a constant frequency 
at 2A = 0.6. Notice also the approach of the frequencies towards 0, as per 
the discussed collision with the linear limit (eigenfunctions) of the problem. 


also preserved for any 7; the relevant eigenfrequency, which 
persists under variations of 7 , can be discerned in the left panel 
of Fig. [ 2 ] For the rest of the spectrum we note that according 
to ( 25 | , eigenfrequencies with nonzero imaginary part can only 
be born in the interval ( — 1 — |A|, 1 + |A|). Here, however, all 
the eigenfrequencies remains inside this interval for all 7 , as 
illustrated in Fig. [ 2 ] solely tending towards 0 , as 7 approaches 
7pt- 

Notice that the VT transition is caused by the nonlinear 
solutions colliding with (or degenerating into) linear modes. 
This fact can also be confirmed in the charge versus 7 plot 
of Fig. [ 3 j where the charge (norm) and energy tend to zero 
(while the width of the solution diverges) when the transition 
point is reached (actually, we have not been able to reach 
this point exactly, as the soliton width increases drastically 
when approaching this point). It should be noted here that 
this is a distinct phenomenology (again) in comparison to the 
NLS counterpart of the model. In the latter, typically at the 
VT -phase transition a stable (center) and an unstable (saddle) 
solution collide and disappear in a saddle-center bifurcation. 
Here, a fundamentally different scenario arises through the 
degeneration of the nonlinear modes into linear ones. In the 
bottom panel of Fig. [ 3 ] we provide two-parameter diagrams 
of the relevant solutions as a function of the frequency A and 
the gain-loss parameter 7 . The dependences strongly suggest a 
“combined” monoparametric dependence on y 2 +A 2 , although 
we have not been able to analytically identify solutions bearing 
this dependence, except in the limit of m = 0; see below. 
For this reason (the strong similarity of the dependence of 
monoparametric plots on A and 7 ), we do not show separately 
the dependence on A or fixed 7 . 

We now turn to the consideration of the dynamical evolution 
of solitons past the VT transition point. As indicated above, 
given their generic stability for 7 < 7 ^ 7 -, we do not consider 




Y Y 



Y 




Fig. 3. The standing waves’ energy E, width W and charge Q dependence for 
A = 0.3 and A = 0.8 are shown in the first three panels. The bottom panels 
present a two-parameter diagram of the dependence of the energy (bottom 
left) and width on the frequency A and the gain-loss parameter 7 . 


the latter case. In the case of 7 > 7^77 we have firstly taken 
as initial condition the soliton for A = 0.8 and 7 = 70 = 0.59 
in the simulation with 7 = 7 S > 7pt = \/l — A 2 = 0 . 6 . We 
observe that if 7 S is close enough to 7pt (he., for a “shallow” 
quench), the density oscillates with a frequency that decreases 
with 7 S — 7 pt (see Fig. [ 4 ]). Notice that the charge of the 
new soliton is always higher than the charge of the initial one 
(see the oscillations of Fig. [ 4 ]) and that the maximum charge 
increases with 7 S . Interestingly, in all of these case examples 
we find that the (7-independent) energy is very well conserved 
as is shown in Fig [ 5 ] When the maximum charge is above a 
threshold (this occurs for > 0 . 995 , i.e., for a deep quench), 
the frequency of the new soliton tends to zero and the solution 
starts to grow indefinitely as shown in Fig. [6] If a smaller 
value of 70 is taken, the same phenomenology persists, but the 
indefinite growth emerges for a smaller value of 7 s . In Fig. [ 5 ] 
we have confirmed that both the energy conservation law and 
the moment equation (10) for the power are satisfied in our 
dynamics. The same is true for the case of Fig. [6] where the 
charge grows exponentially (in the case shown in the figure, for 
which 7 s = 1, as ~ exp( 0 . 088 t); although the characteristic 
growth rate depends on 7 S ). Here, the soliton does not collapse, 
as its shape and width are preserved during the growth. Again, 
this type of growth appears to be very different than, say, the 
collapse in the Hamiltonian NLS model 0 In the latter, the 
width decreases and the amplitude increases, whereas here the 
entire solution grows without changing its spatial distribution. 
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Let us mention that the waveform with oscillating charge 
is fairly generic when the quench is not sufficiently deep 
to cause an exponential growth. Remarkably, such solitary 
waves with oscillating charge can be attained by performing an 
SU(1,1) transformation to a standing wave soliton and have 
been predicted in [36]]. This type of solution for a standing 
wave of frequency A is, in fact, intrinsically connected to the 
invariance of the frequency 2A in the spectrum jl9]. More 
specifically, these bi-frequency, oscillating charge coherent 
structures [which can be dubbed as SU(1,1) solitons] are of 
the form: 


U ( x , t)=a-u(x ) exp(— iAt) — icy + u* (x) exp(iAt), (12) 
V (x, t)=a-v(x) exp(—iAt) — ia+u* (x) exp (iAt), (13) 


with 


a± e C, |<a_| 2 — |ct+| 2 = 1- 


(14) 


In this case, {u(x),v(x)} is the standing wave solution 
with frequency A. Consequently, the charge oscillates with 
a frequency 2A as long as 7 T 0. There is an SU(1,1) family 
of solutions for each value of 7 and A which fulfills the same 
equations that the standing wave solutions 0 satisfy. As a 
result, when 7 ^ / 70 , an SU(1,1) solution with 7 = 7 S 
is apparently dynamically manifested. Since these periodic 
SU(1,1) solutions and the standing wave solutions only exist 
for 7 < m = 1, it is natural to expect that there are no 
nontrivial fixed points for the dynamics for 7 > 1, hence 
giving rise to the observed growth dynamics. 

We have confirmed that the dynamics observed, e.g., in 
Fig. [4] corresponds to SU(1,1) solutions. For instance, for the 
soliton with 70 = 0.59, A = 0.8, when initializing it for 
7 S = 0.9, it spontaneous ly g ives rise to an oscillatory state of 
the above form of Eq. (12) with A = 0.422, ct_ = 1.0847 
and <a + = 0.4201. On the other hand, using a numerically 
exact (up to a prescribed tolerance) solution of our fixed point 
iteration scheme with a given frequency A (for a desired 
7 ), we can select values of a_ and a+ and the exact form 
of Eq. ff2| ) in order to construct, at will, such bi-frequency 
SU(1,1) solutions. An example of this form is shown in the 
panels of Fig. [ 7 ] (even for 7 = 0) for A = 0.5, ct_ = 1.0500 
and a+ = 0.3202. 


B. Massless NLD equation with k = 1 

Finally, we now turn to the analysis of the solutions in the 
special case of m = 0. This case has the peculiarity that it 
features exact solutions available in analytical form even for 
7^0. Due to the functional form of those solutions, we had 
to use a Chebyshev collocation scheme with inhomogenous 
Dirichlet boundary conditions (to numerically identify their 
stability). As we will see below directly from their functional 
form, these solutions are not VT -symmetric, hence the nonlin¬ 
earity becomes responsible for the breaking of VT -symmetry. 
This, in turn, implies that the linear stability eigenvalues do not 
come in quartets [but rather only in conjugate pairs i.e., in pairs 
with opposite Re(cc)]. Furthermore, although the solutions are 
stable for 7 = 0, they become immediately unstable, as soon 
as 7 7 ^ 0. 



Fig. 4. Top two rows: Dynamics for A = 0.8 using as initial condition 
the soliton with 7 = 0.59, but evolving Eq. {lj for 7 = j s , shown in 
the respective panels. Each panel shows the contour plot of the space-time 
evolution of the soliton density. Bottom two rows: same as the top but with 
each panel showing the time evolution of the soliton charge. The dashed line 
corresponds to the charge of the initial condition. 


There are two main exact solutions for a given value of 7 : 


and 


Ul(x) 

Vl(x) 


^== [(-A + ptanh(px)) 
^-2= [(-A - pttmh(px)) 


17], 
17], 


(15a) 

(15b) 


u 2 {x)= 

V 2 (x)= 


2 V=gA 

1 


[(p + Atanh(px)) + rytanh (px)\ , (16a) 
[(— p + Atanh (px)) + i 7 tanh(px)](J 6 b) 


with p = t/ 7 2 + A 2 . This pair of solutions can be transformed 
to another pair by virtue of the transformation {u,v} -A 
as Eq. Q is invariant under transformations 
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Fig. 5. Dynamics for A = 0.8 using as initial condition the soliton with 
7 o = 0.59. Each panel shows the time evolution of the time derivative of 
the soliton charge and compares it with the right hand side of Eq. (To); the 
bottom panels for each example of 7 illustrate the time dependence of the 
energy fluctuations. 




7 = 0.4 


7 = 0.8 




Fig. 6 . Dynamics for A = 0.8 and 7 S = 1 when the soliton with 70 = 0.59 is 
used as initial condition. The left panel shows both components of the spinors 
at different times. The right panel shows the total charge as a function of time. 
The lower panel compares the time derivative of the charge and the right hand 
side of Eq. (To) and shows the time dependence of the energy fluctuations. 
One can clearly observe the (spatially independent) exponential growth of the 
waveform. 


Fig. 7. Dynamics for SU(1,1) solitons with A = 0.5, ol- = 1.0500 and 
= 0.3202. Each panel shows the time evolution of the soliton density. 
The bottom four panels show the temporal evolution of the charge for each 
of the examples in the top panels. 


{U,V} -A { V *, — U*} when m — 0. This transformation in¬ 
troduces the change uo -A co* into the stability eigenfrequencies 
spectrum. It is worth noticing that both solutions 03 and © 
correspond to the same density, i.e. 

\ui(x)\ 2 + \vi(x)\ 2 = \u 2 (x) I 2 + \v 2 (x)\ 2 (17) 

Figure [8] shows the spectral plane for both solutions of 
the massless equation for 7 = 0.01 and A = —0.9. It can 
be observed that the stability eigenfrequencies for the first 
solution are the complex conjugates of the ones of the second 
solution. Nevertheless, in neither case does the spectrum 
present a four-fold symmetry, given the nonlinearity-induced 
breaking of PT-symmetry. 



Fig. 8 . Spectral plane of the exact analytical solutions of Eqs. for 

m = 0, A = —0.9 and 7 = 0.01. The spectra of the two solutions are 
conjugate, but both reflect their non -VT -symmetric nature in the absence of 
eigenfrequency quartets within the spectral plane. 
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Fig. 9. Dynamics of the soliton fj~5) form = 0, A = —0.9 and 7 = 
0.01. Left panel shows the density plot (darker color means deeper amplitude) 
whereas right panel compares the initial density with that of an intermediate 
time. 


We have performed simulations in order to observe the 
effect of unstable modes for both {ui,vi} and {^ 2 ,^ 2 } 
solutions. Those simulations have been performed using a 
Chebyshev collocation method with Neumann boundary con¬ 
ditions using a 2nd-3rd order Runge-Kutta integrator sup¬ 
plemented by a trapezoidal rule and the backward differen¬ 
tiation formula of order 2 (TR-BDF2 algorithm) (37j. We 
have used as initial condition the standing wave perturbed 
along its unstable eigenmode corresponding to the imaginary 
eigenfrequency direction. As a result of the simulation (see 
Fig. we can observe that the soliton dip starts to slowly 
move along the system with constant speed and preserving its 
shape whereas a counter-propagating pair of small amplitude 
waves is emitted. This behaviour is found for both {ui,vi} 
and {^ 2 ,^ 2 } solutions. If the standing waves were perturbed 
following the direction of the eigenmode corresponding to the 
complex eigenvalue pair, the density dip would remain at rest 
and radiation in form of a staggered perturbation is emitted; 
in fact, the eigenmode is spurious because its spatial profile 
possesses a zigzag tail which is not compatible with solutions 
in the continuum limit, hence we do not consider it further 
here. 


IV. Conclusions and Future Challenges 

In the present work, we introduced a VT -symmetric exten¬ 
sion of the Gross-Neveu (or Soler) model as a prototypical 
example of a PT-symmetric nonlinear Dirac equation (PT- 
NLDE). We have found that the model possesses a num¬ 
ber of remarkable and previously unexplored (some even 
in the standard Hamiltonian case) characteristics. The VT- 
symmetric variant used here involved the Dirac matrix 75 
(whose role in the case of spinors with two components was 
played by the Pauli matrix cu) that, as we argued, introduced 
an unusual form of a VT -symmetric setting where for our 
two-component model, it was not transparent that one field 
bears gain and the other loss, but both, in principle, carry 
both gain and loss. Perhaps this special nature of the problem 
is responsible for the remarkable feature that the Hamiltonian 
form of the energy (which is 7 -independent) is still conserved 
in the presence of the gain-loss parameter 7 (while the power 
is not). On the other hand, the nonlinear solutions identified 
here presented an unusual VT phase transition (for the case 
of unit mass) whereby they degenerated into the linear limit 
of the problem. The solutions (which presented an implicit 
mono-parametric dependence on y 2 + A 2 ) were found to be 


stable everywhere within the regime of exact PT-symmetry 
hence we attempted to consider dynamics beyond this interval. 
There, we found the prototypical formation of time-periodic 
solutions whose bi-frequency character was attributed to the 
SU(1,1) symmetry of the equation. Importantly such periodic 
solutions exist even in the Hamiltonian limit of the model. 
Beyond a critical value of 7 , the solutions were found to 
feature an unusual (spatially independent) growth. Finally, 
exact analytical waveforms could only be identified in the 
special limit of the massless problem with m — 0. These 
solutions had an unexpected characteristic as well in that the 
nonlinearity enabled their breaking of the VT - symmetry for 
all values of (A, 7 ) for which they were found to exist. 

Clearly, this is an unusual and exceptionally rich class 
of models at the interface of the emerging theme of VT - 
symmetric media and the mathematically highly nontrivial 
realm of the nonlinear Dirac equations. Understanding in more 
depth any one of the above features [the energy independence 
on 7 , the exponential, spatially independent growth, the spec¬ 
trum of the massless solutions or that of the SU(1,1) solutions] 
would already enable significant advances in this nascent field 
of PT-NLDE models. Generalizations including variants of 
the model [such as the integrable Thirring model, or the 
physically motivated self-interacting model (i.e., only featuring 
respectively nonlinear terms of the form \U\ 2 U and |U| 2 U 
in Eq. ([!]))] would be particularly worthwhile to consider, as 
would extensions to other settings such as two-dimensional 
ones. Such extensions are presently under consideration and 
relevant results will be reported in future publications. 
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